/*
 * Phi.cpp
 *
 *  Created on: May 24, 2011
 *      Author: marchi
 */

#include "Phi.h"
array3<Complex> * Phi::filter=NULL;
array3<Complex> * Phi::one_R=NULL;

Phi & Phi::operator=(const Grid & x){
	Grid & My=*this;
	My=x;
	return *this;
}
Phi & Phi::operator=(array3<Complex> & x){
	Grid temp=x;
	Grid & My=*this;
	My=x;
	return *this;
}
Phi & Phi::operator=(array3<double> & x){
	Grid temp=x;
	Grid & My=*this;
	My=x;
	return *this;
}
Phi & Phi::operator=(const double &  x){
	Grid temp=x;
	Grid & My=*this;
	My=x;
	return *this;
}


Phi & Phi::operator()(const Rho & myrho){
	array3<Complex> & Mfilter=*filter;
	array3<Complex> & Mone_R=*one_R;

	int nfx,nfy,nfz;
	nfx=(nnx % 2 == 0)? nnx/2: nnx/2+1;
	nfy=(nny % 2 == 0)? nny/2: nny/2+1;
	nfz=(nnz % 2 == 0)? nnz/2: nnz/2+1;

	array3<Complex> ro(nnx,nny,nnz);


	for(unsigned i=0;i<nnx;i++)
		for(unsigned j=0;j<nny;j++)
			for(unsigned k=0;k<nnz;k++){
				ro[i][j][k].real()=myrho[i][j][k];
				ro[i][j][k].imag()=0.0;
			}

	fft3d Forward3(-1,ro);
	fft3d Backward3(1,ro);

	Forward3.fft(ro);

	double mw1,mw2,mw3,mw,fact;
	int ia,ja,ka;
	int nx0=static_cast<int>(nnx);
	int ny0=static_cast<int>(nny);
	int nz0=static_cast<int>(nnz);
	double DV=Rho::getDV();

	for(int i=0;i<nx0;i++){
		ia=(i<nfx)?i : i-nx0;
		for(int j=0;j<ny0;j++){
			ja=(j<nfy)?j : j-ny0;
			for(int k=0;k<nz0;k++){
				ka=(k<nfz)?k : k-nz0;
				mw1=oc[XX][XX]*ia+oc[XX][YY]*ja+oc[XX][ZZ]*ka;
				mw2=oc[YY][XX]*ia+oc[YY][YY]*ja+oc[YY][ZZ]*ka;
				mw3=oc[ZZ][XX]*ia+oc[ZZ][YY]*ja+oc[ZZ][ZZ]*ka;
				mw=mw1*mw1+mw2*mw2+mw3*mw3;
				if(one_R) {
					ro[i][j][k]=Mone_R[i][j][k]*ro[i][j][k]*DV;
				}
				else {
					if(!j && !j && !k) fact=0.0;
					else fact=1.0/(mw*M_PI);
					ro[i][j][k]=fact*ro[i][j][k];
				}
				if(filter) ro[i][j][k]=ro[i][j][k]*Mfilter[i][j][k];
			}
		}
	}

	Backward3.fftNormalized(ro);
	*this=ro;
	return *this;
}

void Phi::MakeFilter(int mx,int my, int mz){
	if(!filter) {
		filter=new array3<Complex>;
		filter->Allocate(nnx,nny,nnz);
	}
	array3<Complex> & Mfilter=*filter;

	std::cout<< Mfilter.Nx() << std::endl;
	fft3d Forward3(-1,Mfilter);

	double fact=1.0/(double) ((2*mx+1)*(2*my+1)*(2*mz+1));

	for(int i=-mx;i<=mx;i++){
		int ia=(i>=0)?i:i+nnx;
		for(int j=-mx;j<=mx;j++){
			int ja=(j>=0)?j:j+nny;
			for(int k=-mx;k<=mx;k++){
				int ka=(k>=0)?k:k+nnz;
				Mfilter(ia,ja,ka)=fact;
			}
		}
	}
	Forward3.fft(Mfilter);
	for(unsigned int i=0;i<nnx;i++)
		for(unsigned int j=0;j<nny;j++)
			for(unsigned int k=0;k<nnz;k++)
				Mfilter[i][j][k].imag()=-Mfilter[i][j][k].imag();
}

void Phi::MakeOne_R(){
	try{
		if(nnx+nny+nnz==0) throw "Should initialize Metrics first";
		}
	catch(const char * s){
		std::cout << s << std::endl;
		exit(1);
	}

	if(!one_R) {
		one_R=new array3<Complex>;
		one_R->Allocate(nnx,nny,nnz);
	}

	array3<Complex> & Mone_R=*one_R;
	int nx0=static_cast<int>(nnx);
	int ny0=static_cast<int>(nny);
	int nz0=static_cast<int>(nnz);

	int nfx=(nx0 % 2 == 0)? nx0/2: nx0/2+1;
	int nfy=(ny0 % 2 == 0)? ny0/2: ny0/2+1;
	int nfz=(nz0 % 2 == 0)? nz0/2: nz0/2+1;

	matrix & CO=co;
	fft3d Forward3(-1,Mone_R);
	double ia,ja,ka,rx,ry,rz,rs;
	rx=CO[XX][XX]*(0.5/static_cast<double>(nx0));
	ry=CO[YY][YY]*(0.5/static_cast<double>(ny0));
	rz=CO[ZZ][ZZ]*(0.5/static_cast<double>(nz0));
	double fact=1.0/sqrt((rx*rx+ry*ry+rz*rz)/3.0);

	for(int i=0;i<nx0;i++){
		ia=((i<nfx)?i : i-nx0)/static_cast<double>(nx0);
		for(int j=0;j<ny0;j++){
			ja=((j<nfy)?j : j-ny0)/static_cast<double>(ny0);
			for(int k=0;k<nz0;k++){
				ka=((k<nfz)?k : k-nz0)/static_cast<double>(nz0);
				rx=CO[0][0]*ia+CO[0][1]*ja+CO[0][2]*ka;
				ry=CO[1][0]*ia+CO[1][1]*ja+CO[1][2]*ka;
				rz=CO[2][0]*ia+CO[2][1]*ja+CO[2][2]*ka;
				rs=sqrt(rx*rx+ry*ry+rz*rz);
				Mone_R[i][j][k].real()=rs? 1.0/rs: fact;
				Mone_R[i][j][k].imag()=0.0;
			}
		}
	}

	Forward3.fft(Mone_R);
}
